/**
 * 
 */
package cn.cityhouse.avm.regression;

import java.util.Locale;

import org.apache.commons.math3.analysis.function.Exp;
import org.apache.commons.math3.analysis.function.Sqrt;
import org.apache.commons.math3.distribution.FDistribution;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.exception.OutOfRangeException;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.MatrixDimensionMismatchException;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.QRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.SingularMatrixException;
import org.apache.commons.math3.linear.SingularValueDecomposition;
import org.apache.commons.math3.stat.StatUtils;
import org.apache.commons.math3.stat.descriptive.moment.SecondMoment;
import org.apache.commons.math3.stat.regression.AbstractMultipleLinearRegression;

import cn.cityhouse.avm.stat.SummaryStats;
import cn.cityhouse.avm.stat.WeightedCovMatrices;
import cn.cityhouse.avm.utils.DataFrame;
import cn.cityhouse.avm.utils.GISUtils;
import cn.cityhouse.avm.utils.MatrixUtil2;
import cn.cityhouse.avm.utils.Out;

/**
 * @author <a href="mailto:wangcl@lreis.ac.cn">ChenLiang Wang</a>
 *
 */
public class AVMGWRegression extends AbstractMultipleLinearRegression implements AVMAbstractRegression {

	public enum enum_Kernel {
		GUASS, BISQUARE, EXP
	}

	/**
	 * residuals array of the regression, the default value is {@value}
	 */
	double[] residuals = null;

	/**
	 * R-squared variable member for the regression, the default value is
	 * {@value}
	 */
	double rsquare = Double.NaN;
	/**
	 * The adjusted r-squared variable member for the regression,the default
	 * value is {@value}
	 */
	double adj_rsquare = Double.NaN;
	/**
	 * The
	 */
	double d_sigma = Double.NaN;
	int int_k = -1;// the counter of paramters
	double[] d_y = null;
	double[][] d_x = null;
	/**
	 * Residual Sum of Squares of GWR,the default value is {@value} calculated
	 * by calculateResidualSumOfSquares()
	 */
	double d_SSR = Double.NaN;

	/**
	 * Total Sum of Squares of GWR,the default value is {@value} calculated by
	 * calculateTotalSumOfSquares()
	 */
	double d_SST = Double.NaN;
	/**
	 * Sum of Squared Residuals of GWR,the default value is {@value} calculate
	 * by d_SST - d_SSR
	 */
	double d_SSE = Double.NaN;
	/**
	 * AIC value of GWR
	 */
	double d_AIC = Double.NaN;
	double d_AICc = Double.NaN;
	double[] d_fitted_y = null;
	double vary = Double.NaN;

	/**
	 * enum type indicating which kernel function is used for GWR
	 */
	enum_Kernel em_kernel = enum_Kernel.GUASS;
	boolean isAdapt = false;
	boolean isLonglat = false;
	double[][] pts_data = null;
	double[][] pts_regs = null;
	String[] coord_name = null;
	GWR_dist gDist = null;
	AVMOLRegression olr = null;
	double d_bw = 1e8;
	RealMatrix rlm_beta = null;
	RealMatrix rlm_beta_se = null;
	RealMatrix rlm_beta_tv = null;
	RealVector rlv_fitted = null;
	RealMatrix rlm_Q = null;
	RealMatrix rlm_S = null;
	RealMatrix rlm_Ci = null;
	RealVector rlv_cv = null;
	RealVector rlv_weight_i = null;
	RealMatrix rlm_distance = null;
	double d_trS = Double.NaN;
	double d_trSts = Double.NaN;
	private String[] x_name = null;
	private String y_name = null;
	private boolean isEstimated = false;

	public AVMGWRegression() {
		// TODO Auto-generated constructor stub
		this.olr = new AVMOLRegression();
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see
	 * cn.cityhouse.avm.regression.AVMRegression2#calculateAdjustedRSquared()
	 */
	@Override
	public double calculateAdjustedRSquared() {
		// TODO Auto-generated method stub
		final double n = getX().getRowDimension();
		if (isNoIntercept()) {
			return 1.0d - (1.0d - calculateRSquared()) * (n / (n - getX().getColumnDimension()));
		} else {
			return 1.0d - (calculateResidualSumOfSquares() * (n - 1.0d))
					/ (calculateTotalSumOfSquares() * (n - getX().getColumnDimension()));
		}
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see
	 * org.apache.commons.math3.stat.regression.AbstractMultipleLinearRegression
	 * #calculateBeta()
	 */
	@Override
	protected RealVector calculateBeta() {
		// TODO Auto-generated method stub
		// RealVector rlv=new
		// ArrayRealVector(this.olr.estimateRegressionParameters());
		return calculateBeta(this.rlv_weight_i);
	}

	private RealVector calculateBeta(RealVector rlv_weight_i) {
		// TODO Auto-generated method stub
		// RealVector rlv=new
		// ArrayRealVector(this.olr.estimateRegressionParameters());
		RealMatrix rlm_x = this.getX();
		RealMatrix rlm_xwi = MatrixUtil2.productVectorByCol(rlm_x, rlv_weight_i);
		RealMatrix rlm_tXwi = rlm_xwi.transpose();
		DecompositionSolver solver = new QRDecomposition(rlm_tXwi.multiply(rlm_x)).getSolver();
		RealMatrix rlm_Ci = null;
		RealVector rlv_betai = null;
		RealVector rlv_y = this.getY();
		rlm_Ci = solver.getInverse();
		rlm_Ci = rlm_Ci.multiply(rlm_tXwi);
		this.rlm_Ci = rlm_Ci;
		rlv_betai = rlm_Ci.operate(rlv_y);
		return rlv_betai;
	}

	/**
	 * @param rlv_weight_i
	 *            gw weight
	 * @return CV value
	 */
	private RealVector calculateBeta_CV(RealVector rlv_weight_i_, int i_CV) {
		// gwr.res$lm$x*www.i
		RealVector rlv_weight_i = rlv_weight_i_.copy();
		rlv_weight_i.setEntry(i_CV, 0);// neglect the current point to compute
										// CV score
		RealVector rlm_betai_cv = calculateBeta(rlv_weight_i);
		return rlm_betai_cv;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see
	 * org.apache.commons.math3.stat.regression.AbstractMultipleLinearRegression
	 * #calculateBetaVariance()
	 */
	@Override
	protected RealMatrix calculateBetaVariance() {
		// TODO Auto-generated method stub
		return null;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#calculateFittedValue()
	 */
	@Override
	public RealVector calculateFittedValue() {
		// TODO Auto-generated method stub
		int nrow = olr.getNrow();
		int int_k = olr.getK();
		this.rlm_beta = new Array2DRowRealMatrix(nrow, int_k);
		this.rlm_beta_se = new Array2DRowRealMatrix(nrow, int_k);
		this.rlm_beta_tv = new Array2DRowRealMatrix(nrow, int_k);
		this.rlm_S = new Array2DRowRealMatrix(nrow, nrow);
		this.rlv_fitted = new ArrayRealVector(nrow);
		double[] d_dist_i = null;
		rlv_weight_i = null;
		rlv_cv = new ArrayRealVector(nrow);
		for (int i = 0; i < nrow; i++) {
			if (this.rlm_distance != null && this.rlm_distance.getRowDimension() > 0) {
				d_dist_i = this.rlm_distance.getRow(i);
			} else {
				gDist.focus = i;
				d_dist_i = gDist.computeDist();
			}
			rlv_weight_i = computeGWweight(d_dist_i, this.em_kernel, d_bw);
			this.rlm_beta.setRowVector(i, calculateBeta());
			// if(i==0) Out.printArr(rlm_beta.getRow(i),"|%40.30f");
			this.rlm_S.setRowMatrix(i, getX().getRowMatrix(i).multiply(rlm_Ci));
			RealMatrix tmp = this.rlm_Ci.multiply(this.rlm_Ci.transpose());
			RealVector tmp_diag = MatrixUtil2.getDiagVector(tmp);
			this.rlm_beta_se.setRowVector(i, tmp_diag);
			// this.rlv_fitted.setEntry(i,
			// this.getX().getRowVector(i).dotProduct(this.rlm_beta.getRowVector(i)));
		} // end of nrow
			// Out.printArr(rlm_beta.getRow(0),"|%40.30f");
		this.d_trS = this.rlm_S.getTrace();
		this.d_trSts = MatrixUtil2.getSumOfSquares(rlm_S);
		this.d_SSR = calculateResidualSumOfSquares();
		this.rlv_fitted = rlm_S.operate(this.getY());
		return this.rlv_fitted;
	}

	public double calculateResidualSumOfSquares() {
		// final RealVector residuals = new ArrayRealVector(this.residuals);
		int nrow = this.getNrow();
		RealMatrix rlm_Q1 = MatrixUtils.createRealIdentityMatrix(nrow).subtract(this.rlm_S);
		this.rlm_Q = rlm_Q1.transpose().multiply(rlm_Q1);
		RealMatrix rlm_y = new Array2DRowRealMatrix(this.getY().toArray());
		return rlm_y.transpose().multiply(rlm_Q).multiply(rlm_y).getTrace();
		// No advertised DME, args are valid
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#calculateRSquared()
	 */
	@Override
	public double calculateRSquared() {
		// TODO Auto-generated method stub
		return 1.0d - calculateResidualSumOfSquares() / calculateTotalSumOfSquares();
	}

	/**
	 * Calculate total sum of squares
	 * 
	 * @return total sum of squares
	 */
	public double calculateTotalSumOfSquares() {
		if (isNoIntercept()) {
			return StatUtils.sumSq(this.getY().toArray());
		} else {
			return new SecondMoment().evaluate(this.getY().toArray());
		}
	}

	/**
	 * Compute geographically-weight
	 * 
	 * @param d_dist_i
	 *            distance array
	 * @param enum_kernel_type
	 *            kernel function type
	 * @param d_bw
	 *            bandwidth
	 * @return weight vector for GWR
	 */
	public RealVector computeGWweight(double[] d_dist_i, enum_Kernel enum_kernel_type, double d_bw) {
		RealVector rlv_dist = new ArrayRealVector(d_dist_i);
		RealVector rlv_weight = new ArrayRealVector(d_dist_i.length);
		if (isAdapt) {
			throw new UnsupportedOperationException(
					"Currently,Adaptive Kernel Weighted Methods are not supported.Sorry for that.");
			// switch (enum_kernel_type) {// need to be implemented in the
			// future
			// case GUASS:
			// //rlv_weight = rlv_dist.ebeMultiply(rlv_dist).mapDivide(-2 * d_bw
			// * d_bw).map(new Exp());
			// break;
			// default:
			// break;
			// }
		} else {// Currently , we can only perform fixed bandwidth weighting
				// function
			switch (enum_kernel_type) {
			case GUASS:
				rlv_weight = rlv_dist.ebeMultiply(rlv_dist).mapDivide(-2 * d_bw * d_bw).map(new Exp());
				break;
			case BISQUARE:
				boolean is_Within = false;
				double vdist_cur = Double.NaN;
				for (int i = 0; i < rlv_dist.toArray().length; i++) {
					vdist_cur = rlv_dist.getEntry(i);
					is_Within = vdist_cur < d_bw;
					if (is_Within) {
						double d_weight = 1.0d - vdist_cur * vdist_cur / d_bw / d_bw;
						rlv_weight.setEntry(i, d_weight * d_weight);
					} else {
						rlv_weight.setEntry(i, 0d);
					}
				}
				break;
			case EXP:
				rlv_weight = rlv_dist.mapDivide(-d_bw).map(new Exp());
				break;
			default:
				break;
			}
		}
		return rlv_weight;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#estimateAIC()
	 */
	@Override
	public double estimateAIC() {
		// TODO Auto-generated method stub
		int nrow = this.getNrow();
		double d_sigma_hat2 = d_SSR / nrow;
		return nrow * Math.log(d_sigma_hat2) + nrow * Math.log(2.0d * Math.PI) + nrow + d_trS;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#estimateAICc()
	 */
	@Override
	public double estimateAICc() {
		// TODO Auto-generated method stub
		int nrow = this.getNrow();
		double d_sigma_hat2 = d_SSR / nrow;
		return nrow * Math.log(d_sigma_hat2) + nrow * Math.log(2.0d * Math.PI)
				+ nrow * ((nrow + d_trS) / (nrow - 2.0d - d_trS));
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#estimateFittedValue()
	 */
	@Override
	public double[] estimateFittedValue() {
		// TODO Auto-generated method stub
		return this.calculateFittedValue().toArray();
	}

	/**
	 * estimation function for GWR
	 */
	public void estimateGWR() {
		int nrow = olr.getNrow();
		calculateFittedValue();
		this.residuals = this.estimateResiduals();
		double d_sigma_hat1 = d_SSR / (nrow - 2 * d_trS + d_trSts);
		// RealVector stud_res=new ArrayRealVector(this.residuals);
		// RealVector q_diag=MatrixUtil2.getDiagVector(rlm_Q);
		for (int i = 0; i < nrow; i++) {
			// stud_res.setEntry(i,this.residuals[i]/Math.sqrt(d_sigma_hat1*q_diag.getEntry(i)));
			rlm_beta_se.setRowVector(i, rlm_beta_se.getRowVector(i).mapMultiply(d_sigma_hat1).map(new Sqrt()));
			rlm_beta_tv.setRowVector(i, rlm_beta.getRowVector(i).ebeDivide(rlm_beta_se.getRowVector(i)));
		}
		double d_sigma_hat2 = d_SSR / nrow;
		this.d_AIC = estimateAIC();
		this.d_AICc = estimateAICc();
		double d_edf = d_trSts + nrow - 2.0d * d_trS;
		double d_enp = 2.0d * d_trS - d_trSts;
		this.d_SST = this.calculateTotalSumOfSquares();
		this.rsquare = calculateRSquared();
		this.adj_rsquare = 1.0d - (1.0d - rsquare) * (nrow - 1.0d) / (d_edf - 1.0d);
		isEstimated = true;
	}

	/**
	 * Estimates CV score for d_bw
	 * 
	 * @param d_bw
	 *            bandwidth of gwr
	 * @return the CV score
	 */
	public double estimateGWR_CVscore(double d_bw) {
		int nrow = olr.getNrow();
		int int_k = olr.getK();
		RealMatrix rlm_beta = new Array2DRowRealMatrix(nrow, int_k);
		RealVector rlv_fitted = new ArrayRealVector(nrow);
		double[] d_dist_i = null;
		RealVector rlv_weight_i = null;
		rlv_cv = new ArrayRealVector(nrow);
		for (int i = 0; i < nrow; i++) {
			if (this.rlm_distance != null && this.rlm_distance.getRowDimension() > 0) {
				d_dist_i = this.rlm_distance.getRow(i);
			} else {
				gDist.focus = i;
				d_dist_i = gDist.computeDist();
			}
			rlv_weight_i = computeGWweight(d_dist_i, this.em_kernel, d_bw);
			try {
				rlm_beta.setRowVector(i, calculateBeta_CV(rlv_weight_i, i));
			} catch (SingularMatrixException e) {
				e.printStackTrace();
				this.rlv_cv.setEntry(i, Double.POSITIVE_INFINITY);
				System.out.println("Fixed bandwidth: " + d_bw + "\tCV score: " + Double.POSITIVE_INFINITY);
				return Double.POSITIVE_INFINITY;
			} catch (OutOfRangeException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
				this.rlv_cv.setEntry(i, Double.POSITIVE_INFINITY);
				System.out.println("Fixed bandwidth: " + d_bw + "\tCV score: " + Double.POSITIVE_INFINITY);
				return Double.POSITIVE_INFINITY;
			} catch (MatrixDimensionMismatchException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
				this.rlv_cv.setEntry(i, Double.POSITIVE_INFINITY);
				System.out.println("Fixed bandwidth: " + d_bw + "\tCV score: " + Double.POSITIVE_INFINITY);
				return Double.POSITIVE_INFINITY;
			}
			rlv_fitted.setEntry(i, this.getX().getRowVector(i).dotProduct(rlm_beta.getRowVector(i)));
			this.rlv_cv.setEntry(i, this.getY().getEntry(i) - rlv_fitted.getEntry(i));
		}
		double d_norm_cv = rlv_cv.getNorm();
		System.out.println("Fixed bandwidth: " + d_bw + "\tCV score: " + d_norm_cv * d_norm_cv);
		return (d_norm_cv * d_norm_cv);
	}

	@Override
	public double[] estimateResiduals() {
		// TODO Auto-generated method stub
		return this.getY().subtract(this.rlv_fitted).toArray();
	}

	/**
	 * @return double indicating the bandwidth parameter of GWR
	 */
	public double getD_bw() {
		return d_bw;
	}

	/**
	 * @return string indicating which kernel type is used for GWR
	 */
	public String getEm_kernel() {
		return em_kernel.toString();
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#getK()
	 */
	@Override
	public int getK() {
		// TODO Auto-generated method stub
		return getX().getColumnDimension();
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#getNrow()
	 */
	@Override
	public int getNrow() {
		// TODO Auto-generated method stub
		return getX().getRowDimension();
	}

	@Override
	public RealMatrix getX() {
		return this.olr.getX();
	}

	/**
	 * @return String array of the names for variables
	 */
	public String[] getX_name() {
		return olr.getX_name();
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see
	 * org.apache.commons.math3.stat.regression.AbstractMultipleLinearRegression
	 * #getY()
	 */
	@Override
	public RealVector getY() {
		return this.olr.getY();
	}

	/**
	 * @param data_frame
	 *            a data object containing the variables in the model
	 * @param x_name
	 *            names of the independent variables
	 * @param y_name
	 *            name of the dependent variable
	 */
	private void loadData(DataFrame data_frame, String[] x_name, String y_name) {
		// TODO Auto-generated method stub

		this.olr.loadData(data_frame, x_name, y_name);
	}

	/**
	 * @param data_frame
	 *            a data object containing the variables in the model
	 * @param x_name
	 *            names of the independent variables
	 * @param y_name
	 *            name of the dependent variable
	 * @param coord_name
	 *            names of the coordinates
	 * @param kernel
	 *            String specifing the kernel type of GWR
	 */
	public void loadData(DataFrame data_frame, String[] x_name, String y_name, String[] coord_name, String kernel) {
		if (kernel != null) {
			setEm_kernel(kernel);
		}
		if (coord_name == null || coord_name[0] == null || coord_name[1] == null || coord_name[0].length() < 1
				|| coord_name[1].length() < 1) {
			throw new IllegalArgumentException("Must Provide coordinate!");
		}
		System.out.println("using kernel " + this.em_kernel.toString());
		loadData(data_frame, x_name, y_name);
		this.pts_regs = data_frame.getColumns(coord_name);
		this.isLonglat = GISUtils.isLonglat(new Array2DRowRealMatrix(pts_regs));
		gDist = new GWR_dist(pts_regs);
		rlm_beta = null;
		rlm_beta_se = null;
		rlm_beta_tv = null;
		rlv_fitted = null;
		rlm_Q = null;
		rlm_S = null;
		rlm_Ci = null;
		rlv_cv = null;
		d_trS = Double.NaN;
		d_trSts = Double.NaN;
		rlm_distance = new Array2DRowRealMatrix(this.getNrow(), this.getNrow());
		double[] d_dist_i = null;
		for (int i = 0; i < this.getNrow(); i++) {
			gDist.focus = i;
			d_dist_i = gDist.computeDist();
			this.rlm_distance.setRow(i, d_dist_i);
		}
	}

	/**
	 * @param d_y
	 *            double array of the dependent variable
	 * @param d_x
	 *            double array of the independent variables
	 * @param x_name
	 *            names of the dependent variable
	 * @param y_name
	 *            name of the independent variable
	 */
	private void loadData(double[] d_y, double[][] d_x, String[] x_name, String y_name) {
		this.olr.loadData(d_y, d_x, x_name, y_name);
	}

	/**
	 * @param d_y
	 *            double array of the dependent variable
	 * @param d_x
	 *            double array of the independent variables
	 * @param x_name
	 *            names of the dependent variable
	 * @param y_name
	 *            name of the independent variabl
	 * @param coord_xy
	 *            double array of the coordinates
	 * @param kernel
	 *            String specifing the kernel type of GWR
	 */
	public void loadData(double[] d_y, double[][] d_x, String[] x_name, String y_name, double[][] coord_xy,
			String kernel) {
		if (kernel != null) {
			setEm_kernel(kernel);
		}
		if (coord_xy == null || coord_xy[0] == null || coord_xy[1] == null || coord_xy.length < 1
				|| coord_xy[0].length < 1) {
			throw new IllegalArgumentException("Must Provide coordinate!");
		}
		System.out.println("using kernel " + this.em_kernel.toString());
		loadData(d_y, d_x, x_name, y_name);
		this.pts_regs = coord_xy;
		this.isLonglat = GISUtils.isLonglat(new Array2DRowRealMatrix(pts_regs));
		gDist = new GWR_dist(pts_regs);
		rlm_beta = null;
		rlm_beta_se = null;
		rlm_beta_tv = null;
		rlv_fitted = null;
		rlm_Q = null;
		rlm_S = null;
		rlm_Ci = null;
		rlv_cv = null;
		d_trS = Double.NaN;
		d_trSts = Double.NaN;
		rlm_distance = new Array2DRowRealMatrix(this.getNrow(), this.getNrow());
		double[] d_dist_i = null;
		for (int i = 0; i < this.getNrow(); i++) {
			gDist.focus = i;
			d_dist_i = gDist.computeDist();
			this.rlm_distance.setRow(i, d_dist_i);
		}
	}

	/**
	 * Searching the optimized bandwidth
	 * 
	 * @return the optimized bandwidth parameter
	 */
	public double optimBandwidth() {

		double[][] bbox = GISUtils.getBbox(this.pts_regs);
		double upper = Math.sqrt(Math.pow((bbox[0][1] - bbox[0][0]), 2) + Math.pow((bbox[1][1] - bbox[1][0]), 2));
		double lower = upper / 5000d;
		return searchGoldenSection(lower, upper);
	}

	/**
	 * Predicts coefficients for prediction locations
	 * 
	 * @param rp_locat
	 *            double array of the coordinates for prediction locations
	 * @return RealMatrix indicating the estimated coefficients for the
	 *         prediction locations
	 */
	public RealMatrix predictCoef(double[][] rp_locat) {
		GWR_dist gDist2 = new GWR_dist(this.gDist.getDp_locat(), rp_locat);
		RealMatrix rlm_rp = new Array2DRowRealMatrix(rp_locat);
		RealMatrix rlm_dp = new Array2DRowRealMatrix(this.gDist.getDp_locat());
		int nrow = rlm_rp.getRowDimension();
		System.out.println("Estimating Coefficients at " + nrow + " positions...");
		RealMatrix rlm_dis = new Array2DRowRealMatrix(nrow, rlm_dp.getRowDimension());
		int int_k = getK();
		RealMatrix rlm_beta = new Array2DRowRealMatrix(nrow, int_k);
		double[] d_dist_i = null;
		RealVector rlv_weight_i = null;
		for (int i = 0; i < nrow; i++) {
			gDist2.focus = i;
			d_dist_i = gDist2.computeDist();
			rlm_dis.setRow(i, d_dist_i);
			rlv_weight_i = computeGWweight(d_dist_i, this.em_kernel, d_bw);
			rlm_beta.setRowVector(i, calculateBeta(rlv_weight_i));
		} // end of nrow
			// System.out.println(this.rlv_weight_i);
			// System.out.println(rlv_weight_i);
			// System.out.println(this.rlm_beta);
			// System.out.println(rlm_beta);
		return rlm_beta;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#printSummary()
	 */
	@Override
	public void printSummary() {
		// TODO Auto-generated method stub
		System.out.println("\n----------Results of Global Regression ------------");
		this.olr.printSummary();
		System.out.println("\n---------Results of Geographically Weighted Regression--------");
		System.out.println("\n---------Model calibration information--------");
		System.out.println("Kernel function: " + this.getEm_kernel());
		System.out.println(this.isAdapt ? "Adaptive " : "Fixed " + " bandwidth: " + this.d_bw);
		System.out.println("------------Summary of GWR coefficient estimates:-----------");

		for (int i = 0; i < this.getK(); i++) {
			System.out.println(this.getX_name()[i]);
			SummaryStats.summaryStats(this.rlm_beta.getColumnVector(i).toArray());
		}
		System.out.println("Residuals");
		SummaryStats.summaryStats(this.residuals);
		System.out.println("------------------Diagnostic information--------------");
		System.out.println("Number of data points:\t" + this.olr.getNrow());
		System.out.println("Effective number of parameters (2trace(S) - trace(S'S)):\t" + (2.0d * d_trS - d_trSts));
		System.out.println("Effective degrees of freedom (n-2trace(S) + trace(S'S)):\t"
				+ (d_trSts + this.getNrow() - 2.0d * d_trS));
		// System.out.println("AICc (GWR book, Fotheringham, et al. 2002, p. 61,
		// eq 2.33):\t"+this.estimateAIC());
		// System.out.println("AIC (GWR book, Fotheringham, et al. 2002,GWR p.
		// 96, eq. 4.22):\t"+this.estimateAICc());
		System.out.printf("AIC: %.6f , AICc: %.6f , BIC: %.6f\n", 
				this.estimateAIC(), this.estimateAICc(),this.estimateBIC());
		System.out.printf("Residual sum of squares:\t%15.6f\n", this.d_SSR);
		System.out.println("Multiple R-squared: " + this.rsquare + "\tAdjusted R-squared: " + this.adj_rsquare);
	}

	/**
	 * <a href="https://en.wikipedia.org/wiki/Golden_section_search">Golden
	 * Section Search</a> algorithm for finding the minimum CV score of
	 * bandwidth *
	 * 
	 * @param lower
	 *            lower bound for searching
	 * @param upper
	 *            upper bound for searching
	 * @return the minimum value
	 */
	public double searchGoldenSection(double lower, double upper) {
		double d_eps = 1e-4;// working value ,can be changed
		double R = (Math.sqrt(5) - 1.0d) / 2.0d; // R = 0.61803398....
		int iter = 0;
		double d_d = R * (upper - lower);
		double xL, xU;
		xL = lower;
		xU = upper;
		double x1, x2;
		if (this.isAdapt) {
			x1 = Math.floor(xL + d_d);
			x2 = Math.round(xU - d_d);
			throw new UnsupportedOperationException();
		} else {
			x1 = xL + d_d;
			x2 = xU - d_d;
		}
		double f1 = estimateGWR_CVscore(x1);
		double f2 = estimateGWR_CVscore(x2);
		double d1 = f2 - f1;
		// ****************** Establish initial value of
		// xopt:*******************
		double xopt;
		if (f1 < f2) {
			xopt = x1;
		} else {
			xopt = x2;
		}
		// *************************** main loop********************************
		int ea = 100;

		while ((Math.abs(d_d) > d_eps) && (Math.abs(d1) > d_eps)) {
			d_d = R * d_d;
			if (f1 < f2) {
				xL = x2;
				x2 = x1;
				if (this.isAdapt) {
					x1 = Math.round(xL + d_d);
				} else {
					x1 = xL + d_d;
				}
				f2 = f1;
				f1 = estimateGWR_CVscore(x1);
			} else {
				xU = x1;
				x1 = x2;
				if (this.isAdapt) {
					x2 = Math.floor(xU - d_d);
				} else {
					x2 = xU - d_d;
				}
				f1 = f2;
				f2 = estimateGWR_CVscore(x2);
			}
			iter++;
			// ################# Establish value of xopt after iteration:
			if (f1 < f2) {
				xopt = x1;
			} else {
				xopt = x2;
			}
			d1 = f2 - f1;
		}
		return xopt;
	}

	/**
	 * Set the bandwidth for GWR
	 * 
	 * @param d_bw
	 *            double value indicating the bandwidth parameter for GWR
	 */
	public void setD_bw(double d_bw) {
		this.d_bw = d_bw;
	}

	/**
	 * Set the kernel function for GWR
	 * 
	 * @param em_kernel
	 *            String specifing the kernel type in GWR
	 */
	public void setEm_kernel(String em_kernel) {
		this.em_kernel = enum_Kernel.valueOf(em_kernel.trim().toUpperCase(Locale.ENGLISH));
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see
	 * org.apache.commons.math3.stat.regression.AbstractMultipleLinearRegression
	 * #setNoIntercept(boolean)
	 */
	@Override
	public void setNoIntercept(boolean noIntercept) {
		// TODO Auto-generated method stub
		if (this.olr == null) {
			throw new NullArgumentException();
		}
		super.setNoIntercept(noIntercept);
		this.olr.setNoIntercept(noIntercept);
	}

	/**
	 * Local collinearity diagnostics for basic GWR
	 * 
	 * @return Local condition numbers
	 */
	public RealVector estimateCondNum() {
		RealVector rlv_local_cn = new ArrayRealVector(this.getNrow());
		int nrow = olr.getNrow();
		int int_k = olr.getK();
		double[] d_dist_i = null;
		rlv_weight_i = null;
		rlv_cv = new ArrayRealVector(nrow);
		RealMatrix corr_mat = new Array2DRowRealMatrix(nrow, (int_k - 1) * int_k / 2);
		RealMatrix vdp_idx = new Array2DRowRealMatrix(nrow, int_k);
		for (int i = 0; i < nrow; i++) {
			if (this.rlm_distance != null && this.rlm_distance.getRowDimension() > 0) {
				d_dist_i = this.rlm_distance.getRow(i);
			} else {
				gDist.focus = i;
				d_dist_i = gDist.computeDist();
			}
			rlv_weight_i = computeGWweight(d_dist_i, this.em_kernel, d_bw);
			double sum_wi = MatrixUtil2.getSum(rlv_weight_i);
			rlv_weight_i.mapDivideToSelf(sum_wi);
			int tag = 0;
			for (int j = 0; j < int_k - 1; j++) {
				for (int k = j + 1; k < int_k; k++) {
					tag++;
					RealMatrix rlm_combine = new Array2DRowRealMatrix(nrow, 2);
					rlm_combine.setColumnVector(0, this.getX().getColumnVector(j));
					rlm_combine.setColumnVector(1, this.getX().getColumnVector(k));
					// wrong
					corr_mat.setEntry(i, tag - 1,
							WeightedCovMatrices.estimateWeightedCov(rlm_combine, rlv_weight_i).getEntry(0, 1));
				}
			}
			RealMatrix xw = new Array2DRowRealMatrix(this.getX().getRowDimension(), this.getX().getColumnDimension());
			for (int r = 0; r < xw.getRowDimension(); r++) {
				for (int c = 0; c < xw.getColumnDimension(); c++) {
					xw.setEntry(r, c, this.getX().getEntry(r, c) * rlv_weight_i.getEntry(r));
				}
			}
			RealVector rlv_colsum_xw = MatrixUtil2.getColSum(MatrixUtil2.getElementWiseProduct(xw, xw))
					.mapToSelf(new Sqrt());
			RealMatrix sdv_x = new Array2DRowRealMatrix(this.getX().getRowDimension(),
					this.getX().getColumnDimension());
			for (int c = 0; c < xw.getColumnDimension(); c++) {
				for (int r = 0; r < xw.getRowDimension(); r++) {
					sdv_x.setEntry(r, c, xw.getEntry(r, c) / rlv_colsum_xw.getEntry(c));
				}
			}
			SingularValueDecomposition svd = new SingularValueDecomposition(sdv_x);
			// RealMatrix svd_u = svd.getU();
			RealVector svd_d = MatrixUtil2.getDiagVector(svd.getS());
			// RealMatrix svd_v = svd.getV();
			vdp_idx.setRowVector(i, MatrixUtil2.getReciprocal(svd_d, svd_d.getEntry(0)));
		} // end of nrow
		rlv_local_cn = vdp_idx.getColumnVector(int_k - 1);
		return rlv_local_cn;
	}
	/**
	 * F tests for a cablibrated GWR model without spatial-nonstationary test
	 */
	public void testF() {
		double V1 = this.d_trS;
		double V2 = this.d_trSts;
		int nrow = this.getNrow();
		double RSSg = this.d_SSR;
		double RSSo = this.olr.getD_SSR();
		double DFo = this.olr.d_df_coef_t;
		double delta1 = nrow - 2 * V1 + V2;
		RealMatrix Q = this.rlm_Q;
		double sigma2_delta1 = RSSg / delta1;
		double sigma2 = RSSg / nrow;
		double odelta2 = MatrixUtil2.getSumOfSquares(MatrixUtil2.getDiagVector(Q));
		double delta2 = StatUtils.sum(MatrixUtil2.getDiagVector(Q.multiply(Q)).toArray());
		double L_delta1 = Q.getTrace();
		double L_delta2 = delta2;
		/////////////// F1 test////////////according to Leung et al. 2000
		double F1 = (RSSg / L_delta1) / (RSSo / DFo);
		double[] F1_DF = new double[] { L_delta1 * L_delta1 / L_delta2, DFo };
		double F1_pv = 0;
		FDistribution f = new FDistribution(F1_DF[0], F1_DF[1]);
		F1_pv = f.cumulativeProbability(F1);// lower tail
		System.out.println("---------------F test results of GWR calibration---------------");
		System.out
				.println("---Test for goodness of fit ---------F1&F2 test (Leung et al. 2000),F4 test (GWR book p92)");
		String[] F_name = new String[] { "Numerator DF\t", "Numerator DF\t", "Denominator DF\t", "Pr(>)\t" };
		System.out.printf("%15s", "F statistic:\t");
		// System.out.printf("%15s", "Numerator DF:\t");
		// System.out.printf("%15s", "Denominator DF:\t");
		// System.out.printf("%15s", "Pr(>):\t");
		Out.printArr(F_name);
		System.out.println();
		System.out.printf("%15s", "F1\t");
		Out.printArr(new double[] { F1, F1_DF[0], F1_DF[1] });
		System.out.printf("|%15.6e", F1_pv);
		////////////////// F2 test
		double F2 = ((RSSo - RSSg) / (DFo - L_delta1)) / (RSSo / DFo);
		double[] F2_DF = new double[] { (DFo - L_delta1) * (DFo - L_delta1) / (DFo - 2 * L_delta1 + L_delta2), DFo };
		f = new FDistribution(F2_DF[0], F2_DF[1]);
		double F2_pv = 1 - f.cumulativeProbability(F2);// uppper tail
		System.out.printf("\n%15s", "F2\t");
		Out.printArr(new double[] { F2, F2_DF[0], F2_DF[1] });
		System.out.printf("|%15.6e", F2_pv);
		///////////////////////// F4, the test statistic adopted from the GWR
		///////////////////////// book p92
		double F4 = RSSg / RSSo;
		double[] F4_DF = new double[] { delta1, DFo };
		f = new FDistribution(F4_DF[0], F4_DF[1]);
		double F4_pv = f.cumulativeProbability(F4);// lower tail
		// System.out.println("---F4 test (GWR book p92)");
		// System.out.printf("F4 statistic: %20.10f ",F4);
		// System.out.print("Numerator DF: "+F4_DF[0]+"\tDenominator DF:
		// "+F4_DF[1]);
		// System.out.printf(" Pr(>) :%20.10e\n",F4_pv);
		System.out.printf("\n%15s", "F4\t");
		Out.printArr(new double[] { F4, F4_DF[0], F4_DF[1] });
		System.out.printf("|%15.6e", F4_pv);		
	}
	/**
	 * F tests for a cablibrated GWR model
	 */
	public void testF1234() {
		double V1 = this.d_trS;
		double V2 = this.d_trSts;
		int nrow = this.getNrow();
		double RSSg = this.d_SSR;
		double RSSo = this.olr.getD_SSR();
		double DFo = this.olr.d_df_coef_t;
		double delta1 = nrow - 2 * V1 + V2;
		RealMatrix Q = this.rlm_Q;
		double sigma2_delta1 = RSSg / delta1;
		double sigma2 = RSSg / nrow;
		double odelta2 = MatrixUtil2.getSumOfSquares(MatrixUtil2.getDiagVector(Q));
		double delta2 = StatUtils.sum(MatrixUtil2.getDiagVector(Q.multiply(Q)).toArray());
		double L_delta1 = Q.getTrace();
		double L_delta2 = delta2;
		/////////////// F1 test////////////according to Leung et al. 2000
		double F1 = (RSSg / L_delta1) / (RSSo / DFo);
		double[] F1_DF = new double[] { L_delta1 * L_delta1 / L_delta2, DFo };
		double F1_pv = 0;
		FDistribution f = new FDistribution(F1_DF[0], F1_DF[1]);
		F1_pv = f.cumulativeProbability(F1);// lower tail
		System.out.println("---------------F test results of GWR calibration---------------");
		System.out
				.println("---Test for goodness of fit ---------F1&F2 test (Leung et al. 2000),F4 test (GWR book p92)");
		String[] F_name = new String[] { "Numerator DF\t", "Numerator DF\t", "Denominator DF\t", "Pr(>)\t" };
		System.out.printf("%15s", "F statistic:\t");
		// System.out.printf("%15s", "Numerator DF:\t");
		// System.out.printf("%15s", "Denominator DF:\t");
		// System.out.printf("%15s", "Pr(>):\t");
		Out.printArr(F_name);
		System.out.println();
		System.out.printf("%15s", "F1\t");
		Out.printArr(new double[] { F1, F1_DF[0], F1_DF[1] });
		System.out.printf("|%15.6e", F1_pv);
		////////////////// F2 test
		double F2 = ((RSSo - RSSg) / (DFo - L_delta1)) / (RSSo / DFo);
		double[] F2_DF = new double[] { (DFo - L_delta1) * (DFo - L_delta1) / (DFo - 2 * L_delta1 + L_delta2), DFo };
		f = new FDistribution(F2_DF[0], F2_DF[1]);
		double F2_pv = 1 - f.cumulativeProbability(F2);// uppper tail
		System.out.printf("\n%15s", "F2\t");
		Out.printArr(new double[] { F2, F2_DF[0], F2_DF[1] });
		System.out.printf("|%15.6e", F2_pv);
		///////////////////////// F4, the test statistic adopted from the GWR
		///////////////////////// book p92
		double F4 = RSSg / RSSo;
		double[] F4_DF = new double[] { delta1, DFo };
		f = new FDistribution(F4_DF[0], F4_DF[1]);
		double F4_pv = f.cumulativeProbability(F4);// lower tail
		// System.out.println("---F4 test (GWR book p92)");
		// System.out.printf("F4 statistic: %20.10f ",F4);
		// System.out.print("Numerator DF: "+F4_DF[0]+"\tDenominator DF:
		// "+F4_DF[1]);
		// System.out.printf(" Pr(>) :%20.10e\n",F4_pv);
		System.out.printf("\n%15s", "F4\t");
		Out.printArr(new double[] { F4, F4_DF[0], F4_DF[1] });
		System.out.printf("|%15.6e", F4_pv);
		///////////////////////////// F3 test
		int var_n = this.getK();
		RealMatrix ek = MatrixUtils.createRealIdentityMatrix(var_n);
		RealMatrix iden = MatrixUtils.createRealIdentityMatrix(nrow);
		RealMatrix J = MatrixUtil2.getRealMatrixWithFilledValue(nrow, nrow, 1);
		RealVector Vk2 = new ArrayRealVector(var_n);
		for (int i = 0; i < var_n; i++) {
			RealMatrix rlm_betaiT = this.rlm_beta.getColumnMatrix(i).transpose();
			RealMatrix rlm_iden_J = iden.subtract(J.scalarMultiply(1.0d / nrow));
			RealMatrix rlm_beta_ = rlm_betaiT.multiply(rlm_iden_J).multiply(this.rlm_beta.getColumnMatrix(i));
			Vk2.setEntry(i, rlm_beta_.getTrace() / nrow);
		}
		RealVector gamma1 = new ArrayRealVector(var_n);
		RealVector gamma2 = new ArrayRealVector(var_n);
		RealVector numdf = new ArrayRealVector(var_n);
		RealVector F3 = new ArrayRealVector(var_n);
		RealVector F3_pv = new ArrayRealVector(var_n);
		RealMatrix F3_DF = new Array2DRowRealMatrix(var_n, 2);
		for (int i = 0; i < var_n; i++) {
			RealMatrix B = new Array2DRowRealMatrix(nrow, nrow);
			for (int j = 0; j < nrow; j++) {
				gDist.focus = j;
				// double[] d_dist_j = gDist.computeDist();
				double[] d_dist_j = null;
				if (this.rlm_distance != null && this.rlm_distance.getRowDimension() > 0) {
					d_dist_j = this.rlm_distance.getRow(j);
				} else {
					gDist.focus = j;
					d_dist_j = gDist.computeDist();
				}
				RealVector rlv_weight_j = computeGWweight(d_dist_j, this.em_kernel, d_bw);
				// RealMatrix diag_Wj=new
				// DiagonalMatrix(rlv_weight_j.toArray());
				// RealMatrix
				// rlm_tXwj=this.getX().transpose().multiply(diag_Wj);
				// RealMatrix
				// rlm_hat=MatrixUtils.inverse(rlm_tXwj.multiply(this.getX()));
				// B.setRowMatrix(j,
				// ek.getRowMatrix(i).multiply(rlm_hat).multiply(rlm_tXwj));
				RealMatrix rlm_x = this.getX();
				RealMatrix rlm_xwi = MatrixUtil2.productVectorByCol(rlm_x, rlv_weight_j);
				RealMatrix rlm_tXwi = rlm_xwi.transpose();
				DecompositionSolver solver = new QRDecomposition(rlm_tXwi.multiply(rlm_x)).getSolver();
				RealMatrix rlm_Ci = solver.getInverse().multiply(rlm_tXwi);
				B.setRowMatrix(j, ek.getRowMatrix(i).multiply(rlm_Ci));
			}
			RealMatrix rlm_iden_J = iden.subtract(J.scalarMultiply(1.0d / nrow));
			RealMatrix BJ = B.transpose().multiply(rlm_iden_J).multiply(B).scalarMultiply(1.0d / nrow);
			RealVector rlv_diag_BJ = MatrixUtil2.getDiagVector(BJ);
			gamma1.setEntry(i, StatUtils.sum(rlv_diag_BJ.toArray()));
			gamma2.setEntry(i, MatrixUtil2.getSumOfSquares(rlv_diag_BJ));
			numdf.setEntry(i, gamma1.getEntry(i) * gamma1.getEntry(i) / gamma2.getEntry(i));
			F3.setEntry(i, Vk2.getEntry(i) / gamma1.getEntry(i) / sigma2);
			f = new FDistribution(numdf.getEntry(i), F1_DF[0]);
			F3_pv.setEntry(i, 1 - f.cumulativeProbability(F3.getEntry(i)));// uppper
																			// tail
			F3_DF.setRow(i, new double[] { numdf.getEntry(i), F1_DF[0] });
		}
		System.out.println("\n---Test for variation of each set of parameters(Leung et al. 2000)");
		System.out.printf("%15s", "Coefficients:\t");
		Out.printArr(getX_name());
		System.out.printf("\n%15s", "F3 statistic:\t");
		Out.printArr(F3.toArray());
		System.out.printf("\n%15s", "Numerator DF:\t");
		Out.printArr(F3_DF.getColumn(0));
		System.out.printf("\n%15s", "Denominator DF:\t");
		Out.printArr(F3_DF.getColumn(1));
		System.out.printf("\n%15s", "Pr(>):\t");
		Out.printArr(F3_pv.toArray(), "|%15.6e");
		System.out.println();
	}

	@Override
	public double estimateBIC() {
		// TODO Auto-generated method stub
		int nrow = this.getNrow();		
		double BIC = this.estimateAIC()-2.0d*(this.getK()+1.0d)+
				Math.log(nrow) * (this.getK() + 1.0d);
		return BIC;
	}
}
